Wyzwanie konkursowe polega na predykcji wskaźnika bezrobocia na lata 2022-2023, na podstawie danych z 2000-2021. Analizę będziemy przeprowadzać przy użyciu języka R w programie RStudio, do obliczeń i przewidywań użyjemy m.in. z bibliotek takich jak forecast, prophet, WSTAW POTEM, do wizualizacji - ggplot2, a do utworzenia ostatecznego raportu posłużymy się oprogramowaniem RMarkdown.
Wejściowe dane przyjmują formę szeregu czasowego - realizacji procesu stochastycznego, którego dziedziną jest czas - w tym przypadku po 12 miesięcy z 22 lat. Procesem stochastycznym \((X_t)_{t \in T}\) nazywamy rodzinę zmiennych losowych z pewnej przestrzeni probabilistycznej, przyjmującą wartości z przestrzeni mierzalnej.
Dane wskaźnika bezrobocia w latach 2000-2021 przedstawiają się w następujący sposób:
Z podziałem na lata:
Od razu zauważamy, że dane podlegają dostrzegalnemu, nieliniowemu trendowi - wartość wskaźnika spada na początku roku, w okolicy połowy roku wzrasta, by potem delikatnie spadać, lub utrzymywać się aż do października, a na koniec roku znowu wzrasta. Występuje oczywiście sezonowość o okresie 12 miesięcy
By sprawdzić podejrzenia wynikające z wizualnej obserwacji powyższych wykresów, przeprowadzamy dekompozycję STL. Ma ona wskazać składniki szeregu czasowego - jego trend, sezonowość oraz reszty.
library(forecast)
ts_stl <- ts(df$Wskaznik, frequency = 12, start = c(2000,12))
autoplot(mstl(ts_stl),col=ts_stl)
Powyższa metoda opiera się na przedstawieniu punktów szeregu czasowego (\(y_i, i \in T\)) jako suma komponentów sezonowości \(s_i\), trendu \(t_i\) i reszty \(r_i\): \[y_i = s_i + t_i + r_i\] oraz estymacja owych komponentów. (CITATION NEEDED)
Obserwując surowe dane, widzimy pewną anomalię w roku 2020 - wyraźny skok wartości wskaźnika bezrobocia w kwietniu w wyniku wybuchu pandemii COVID-19. To zaburzenie w danych w większości przypadków obniży jakość prognozy, ponieważ trend w tym okresie zostaje naruszony. Z tym problemem możemy poradzić sobie na kilka sposobów. Istnieje opcja podstawienia średniej wartości wskaźnika z każdego miesiąca do odpowiednich miesięcy z 2020. Można wziąć średnie globalne, lub jedynie z kilku ostatnich lat. Nie ma również większych przeszkód, by zupełnie pominąć ten rok w obliczeniach. Zbadamy także ideę, by użyć danych z lat 2000-2019, by “przewidzieć” wartości z 2020, a następnie dokonywać obliczeń przy użyciu nowych danych z 2020 do predykcji 2022-2023. Dokonamy analizy tych metod w rozdziale 3.
Kolejną metodą przewidywania szeregu czasowego jest Prophet, przedstawiony przez Facebooka.[*] Na wyjściowym szeregu dokonujemy dekompozycji w następujący sposób: \[y(t) = g(t) + s(t) + h(t) + e_t\] W tym przypadku \(g(t)\) jest funkcją trendu reprezentującą nieokresowe zmiany wartości szeregu czasowego, \(s(t)\) jest funkcją zmian okresowych (np. miesięcznych), a \(h(t)\) to efekt świąt, występujących w nieregularnych odstępach czasu. Zakładamy, że błąd \(e_t\) ma rozkład normalny. funkcje \(g(t)\), \(s(t)\), \(h(t)\) są estymowane przy użyciu m. in. logistycznego modelu wzrostu oraz szeregów Fouriera, ściślej opisane w [*].
Podstawiając nasze dane uzyskujemy następujące wyniki predykcji:
library(prophet)
df_prophet <- df
df_prophet$date <- as.Date(paste(df$Rok, df$M.c, "01", sep = "-"))
df_prophet <- subset(df_prophet, select = c(date, Wskaznik))
colnames(df_prophet) <- c('ds', 'y')
model_prophet <- prophet(df_prophet)
future <- make_future_dataframe(model_prophet,
periods = 24,
freq = 'month')
forecast <- predict(model_prophet, future)
forecast[c('ds', 'yhat')] %>%
tail(n = 24)
## ds yhat
## 265 2022-01-01 105.59152
## 266 2022-02-01 100.91219
## 267 2022-03-01 98.81112
## 268 2022-04-01 97.07026
## 269 2022-05-01 97.14200
## 270 2022-06-01 97.91921
## 271 2022-07-01 99.66101
## 272 2022-08-01 99.95366
## 273 2022-09-01 99.73869
## 274 2022-10-01 99.50115
## 275 2022-11-01 101.79561
## 276 2022-12-01 103.00231
## 277 2023-01-01 105.79901
## 278 2023-02-01 101.15992
## 279 2023-03-01 99.35893
## 280 2023-04-01 97.09987
## 281 2023-05-01 97.20263
## 282 2023-06-01 98.16253
## 283 2023-07-01 100.04565
## 284 2023-08-01 100.28066
## 285 2023-09-01 99.94061
## 286 2023-10-01 99.66011
## 287 2023-11-01 101.96982
## 288 2023-12-01 103.25350
prophet_plot_components(model_prophet, forecast)
plot(model_prophet, forecast)
Na pierwszym wykresie widzimy krzywą trendu wyznaczoną przez model, łącznie z przewidzianymi ostatnimi latami wraz z przedziałem ufności (niebieskie pole na końcu krzywej), na drugim rysunku - uśrednione, ogólne roczne zmiany wskaźnika. Ostatni wykres przedstawia dodatkowo porównanie rzeczywistych danych (czarne punkty) z tymi estymowanymi przez Prophet (niebieska linia, wraz z przedziałem ufności).
Zbadajmy teraz jakość predykcji, kiedy pominiemy rok 2020 w rozważaniach:
forecast2[c('ds', 'yhat')] %>%
tail(n = 24)
## ds yhat
## 253 2022-01-01 104.48811
## 254 2022-02-01 99.89273
## 255 2022-03-01 97.55734
## 256 2022-04-01 95.63517
## 257 2022-05-01 95.70601
## 258 2022-06-01 96.53879
## 259 2022-07-01 98.34170
## 260 2022-08-01 98.68818
## 261 2022-09-01 98.51341
## 262 2022-10-01 98.31010
## 263 2022-11-01 100.63690
## 264 2022-12-01 101.82312
## 265 2023-01-01 104.36696
## 266 2023-02-01 99.81617
## 267 2023-03-01 97.72722
## 268 2023-04-01 95.69684
## 269 2023-05-01 95.67828
## 270 2023-06-01 96.49088
## 271 2023-07-01 98.35158
## 272 2023-08-01 98.70011
## 273 2023-09-01 98.48975
## 274 2023-10-01 98.33266
## 275 2023-11-01 100.65959
## 276 2023-12-01 101.92684
prophet_plot_components(model_prophet2, forecast2)
Obserwujemy, że po wyrzuceniu 2020, linia trendu zdecydowanie się wypłaszcza na koniec badanego okresu. Przez to, że model bierze pod uwagę współczynnik świąt, czyli nieregularnych skoków badanej zmiennej, w tym przypadku warto zostawić dane z 2020, zatem ostatecznie bierzemy pierwotną wersję prognozy.
PROPHET <- forecast['yhat'] %>%
tail(n = 24)
rownames(PROPHET) <- 1:24
colnames(PROPHET) <- 'Prophet'
## Prophet
## 1 105.59152
## 2 100.91219
## 3 98.81112
## 4 97.07026
## 5 97.14200
## 6 97.91921
## 7 99.66101
## 8 99.95366
## 9 99.73869
## 10 99.50115
## 11 101.79561
## 12 103.00231
## 13 105.79901
## 14 101.15992
## 15 99.35893
## 16 97.09987
## 17 97.20263
## 18 98.16253
## 19 100.04565
## 20 100.28066
## 21 99.94061
## 22 99.66011
## 23 101.96982
## 24 103.25350
OPIS
TEORETYCZNY
TEJ
METODY
BIBLIOGRAFIA[**]
W celu ewaluacji jakości prognozy, będziemy porównywać ostatnie dane 2 lata z przewidzianymi wartościami na następne 2 lata, korzystając z dwóch metryk:
Po podstawieniu naszych danych otrzymujemy następujące wyniki:
library(forecast)
ts_tbats <- msts(df$Wskaznik, seasonal.periods = 12)
model.ts_tbats <- tbats(ts_tbats)
model.ts_tbats.forecast <- forecast(model.ts_tbats, h = 24)
plot(model.ts_tbats.forecast, main = 'Prognoza TBATS',
ylab = 'Wskaznik')
model.ts_tbats.forecast$mean
## Jan Feb Mar Apr May Jun Jul
## 23 102.29247 98.84214 96.22711 95.65658 95.33452 96.74659 97.75088
## 24 104.10502 100.24764 97.32554 96.53332 96.03627 97.31860 98.21517
## Aug Sep Oct Nov Dec
## 23 98.91372 98.21906 98.66880 100.31739 102.23146
## 24 99.29118 98.52024 98.91193 100.51604 102.39415
df.test <- tail(df$Wskaznik, n = 12 * 2)
ts_tbats.predict <- predict(model.ts_tbats.forecast, df.test)
ts_tbats.test_vs_predicted <- data.frame(df.test, model.ts_tbats.forecast$mean)
matplot(ts_tbats.test_vs_predicted, type = 'l', lty = 1:2, col = 1:2, ylab = 'Przypadek nr 1')
MAPE.error(ts_tbats.test_vs_predicted$df.test,
ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 2.085391
RMSE.error(ts_tbats.test_vs_predicted$df.test,
ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 6.855342
W wykresie prognozy, niebieska linia oznacza przewidziane dane, wraz z przedziałami ufności na poziomie istotności \(0.2\) oraz \(0.05\). Przewidziane wartości zostały wypisane w tabeli pod nim.
Na ostatnim wykresie widzimy porównanie wartości z ostatnich dwóch lat (czarna linia) z przewidzianymi wartościami na następne 2 lata (czerwona przerywana linia)[***]. Obserwujemy tutaj zauważalne niedopasowanie krzywych na początku okresu - jest to następstwo wzięcia pod uwagę odstających wartości z roku 2020. Ostatecznie, liczymy wartości błędów MAPE i RMSE - są one zawyżone z wyżej wymienionego powodu.
Aby zmiejszyć wartości błędów, dokonujemy analizy bez 2020:
## MAPE: 1.168697
## RMSE: 3.977294
W przypadku wyrzucenia danych z 2020, zgodnie z oczekiwaniami, widzimy znaczny spadek błędów MAPE i RMSE. Zwróćmy uwagę, że do ewaluacji prognozy wzięliśmy lata 2019 i 2021 zamiast 2020-2021. Sprawdzimy jeszcze predykcję dla przypadków, gdy zamiast wartości z tego roku podstawimy średnią globalną wartość wskaźnika:
## MAPE: 1.009071
## RMSE: 2.643147
Po raz kolejny, wartości błędu zmiejszyły się. Na koniec zbadajmy zachowanie błędów, gdy zamiast 2020 weźmiemy średni wskaźnik z ostatnich 5 lat:
## MAPE: 1.06453
## RMSE: 3.255959
Najmniejsze wartości błędów otrzymujemy w przypadku podstawienia globalnej średniej za wartości z odstającego roku. Decydujemy jednak uwzględnić ostatni przypadek w ostatecznej predykcji, gdyż nie chcemy brać takich wyników, które będą zbyt zbliżone do wartości z ostatnich dwóch lat, ponieważ po raz kolejny, trend byłby zbyt spłaszczony na końcu okresu. Ostatecznie:
TBATS <- model.ts_tbats4.forecast$mean %>% as.data.frame()
colnames(TBATS) <- 'TBATS'
## TBATS
## 1 102.15098
## 2 98.84227
## 3 96.15523
## 4 95.20574
## 5 94.80567
## 6 96.32381
## 7 97.41151
## 8 98.46648
## 9 97.71175
## 10 98.19079
## 11 99.75118
## 12 101.52969
## 13 103.32886
## 14 99.75306
## 15 96.86344
## 16 95.76631
## 17 95.25199
## 18 96.68639
## 19 97.70474
## 20 98.70353
## 21 97.89990
## 22 98.34201
## 23 99.87406
## 24 101.62973
** De Livera A.M., Hyndman R.J., Snyder R.D., Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American Statistical Association 2011, 106, 1513–1527.
*** Jasiewicz-Badowska M. Analiza i prognozowanie szeregów czasowych z uwzględnieniem sezonowości 2016